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Abstract 

A technique  is  presented  for  rapidly  producing 
unstructured  finite  element  meshes  in  support  of  large- 
scale  remote  sensing  simulations.  These  tetrahedral 
meshes  typically  have  more  than  one  million  elements  and 
more  than  250  thousand  nodes,  and  allow  for  arbitrary 
placement  of  objects  into  the  domain.  Open-source  mesh 
generation  packages  are  used  in  conjunction  with  a 
tetrahedra  element  smoothing  operation  to  achieve  the 
desired  final  meshes.  Meshes  can  be  reproduced  in  less 
than  30  minutes  on  a Cray  XT3  architecture. 

1.  Introduction 

As  part  of  an  ongoing  effort  to  improve  automated 
target  recognition  algorithms  for  remote  sensing 
technologies,  a suite  of  closely  coupled  numerical 
simulators  has  been  developed.  These  models  include 
high-resolution  thermal  and  moisture  transport  finite 
element  models,  as  well  as  solar  (ray  casting)  and 
vegetation  models.  The  resulting  tightly  coupled  system 
allows  heat  transport  simulations  to  be  performed  at  a 
very  fine  scale.  This  setup  is  well  suited  for  the  testing  of 
specific  scenes  with  controlled  environmental  conditions, 
in  particular  meteorological  and  time-of-day  conditions, 
which  might  otherwise  be  difficult  and  time  consuming  to 
reproduce  in  the  field.  However,  this  suite  serves  as  a 
complement  to,  rather  than  a replacement  for,  field  and 
laboratory  testing  of  remote  sensors. 

These  numerical  simulations  require  the  rapid  and 
robust  production  of  unstructured  finite  element  meshes 
of  the  shallow  subsurface  for  the  moisture  and  thermal 
codes.  The  meshes  need  to  include  natural  and  man-made 
objects  to  achieve  realism  and  relevancy,  e.g.,  rocks  and 
soda  cans.  This  work  will  focus  on  the  mesh  generation 
process,  which  has  been  achieved  by  taking  advantage  of 
open-source  (black-box)  mesh  generation  software.  The 


tetrahedra  that  are  output  are  then  post-processed  with  a 
mesh  smoothing  technique  to  ensure  quality  elements  in 
the  final  tetrahedral  mesh.  In  addition,  the  process  allows 
for  the  inserted  objects  to  be  arbitrarily  placed,  i.e., 
buried,  flush  with  or  protruding  from  the  ground  surface. 
A simple  mesh  repair  operation  around  the  objects  is 
utilized  to  avoid  poorly  shaped  tetrahedra  in  these 
regions.  When  desired,  subsurface  soil  regions  can  be 
assigned  as  a post-process  step.  This  ability  extends  to 
statistically  generated  soil  distributions  using  site-specific 
information  obtained  from  soil  samples. 

2.  Countermine  Computational  Testbed 

Collectively,  the  numerical  simulators  and  their  pre- 
and  post-processing  tools  are  called  the  Countermine 
Computational  Testbed.  The  purpose  of  the  testbed  is  to 
represent  the  thermal,  hydrological,  and  meteorological 
processes  that  contribute  significantly  to  the  infrared 
response  of  typical  scenes.  The  testbed  must  produce 
visually  realistic  images,  which  requires  spatial  variability 
at  a resolution  that  is  the  same  or  better  than  that 
measurable  by  the  fielded  sensor  (about  2-4  cm). 

2.1.  Scene  Generation 

A typical  scene  includes  soils,  rocks,  plants,  and 
targets  or  other  objects.  The  soil  volume  normally 
extends  0.5  to  1 m into  the  ground  and  includes  a detailed 
representation  of  the  ground  surface.  Surface  texture 
dictates  the  spatial  variability  in  surface  temperatures, 
especially  at  low  sun  angles  and  with  sparse  vegetation. 
Spatial  distribution  of  soil  materials  also  contributes  to 
temperature  variability.  High  variability  tends  to  make 
target  discrimination  more  difficult. 

For  thermal  infrared  scene  modeling,  the  vegetative 
elements  on  the  surface  may  be  represented  by  either  a 
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plane  parallel  abstraction  or  as  individual  three- 
dimensional  geometric  models.  Because  of  the  sensor 
ground  resolution  requirements,  the  plants  are  modeled  as 
individual  geometric  objects  where  the  stems  and  leaves 
are  described  as  triangular  facets.  Geo-typical 
representations  of  the  plants  are  based  on  in  situ 
measurements.  Thermal  conduction  in  plant  materials  is 
relatively  small  and  the  thermal  analysis  of  the  plants  can 
be  computed  independently  of  the  ground  model.  Thus, 
the  plants  are  not  an  integral  part  of  the  finite  element 
mesh  but  are  merely  placed  on  the  surface  mesh. 

2.2.  Terrain  Physics  Models 

Energy  is  introduced  to  the  scene  by  a ray  casting 
model.  The  ray  caster  uses  the  sun  position,  a discrete 
representation  of  the  hemispherical  sky,  and  hourly 
meteorological  data  to  compute  inbound  radiation 
distributed  among  six  spectral  bands.  These  loadings  are 
applied  to  individual  soil,  vegetation,  rock,  and  target 
facets  in  the  scene.  The  ray  caster  includes  shading  of 
facets  by  other  facets  and  permits  partial  transmission  of 
energy  through  plant  facets.  Facets  reflect  a portion  of 
the  energy  within  each  spectral  band.  The  vegetation 
model  accepts  the  energy  input  from  the  ray  caster  and 
computes  leaf  temperatures  based  on  material  properties, 
air  temperature,  wind  speed,  etc.  In  addition  to  plant  facet 
temperatures,  the  vegetation  model  estimates  stomatal 
resistance  and  the  resulting  water  demand  from  its  root 
system. 

The  soil  moisture  and  thermal  migration  model 
(ADH)t8,101  is  a finite  element  model  that  uses  linear  basis 
functions  on  tetrahedra.  The  model  approximates 
Richards  equation  for  partially  saturated  groundwater 
flow  and  heat  conduction  and  convection  with  moisture- 
content-dependent  properties.  ADH  accepts  rainfall  from 
the  meteorological  data  and  inbound  energy  from  the  ray 
caster  as  some  of  its  boundary  conditions.  Heat  flux  at 
the  ground  surface  also  depends  on  sensible  and  latent 
heat  exchange.  Mesh  refinement  and  coarsening  can  be 
triggered  by  multiphysics  error  indicators.  The  mesh 
refinement  scheme  attempts  to  retain  the  original  element 
aspect  ratio  by  splitting  oldest  edges  first.  The  model 
accommodates  parallel  processing  with  dynamic  load 
balancing. 

The  product  of  the  vegetation  and  soil  models  is 
predicted  temperature  at  each  of  the  nodes  for  the  time 
step.  The  ray  caster  uses  these  nodal  temperatures  to 
compute  emitted  energy,  which  is  added  to  the  reflected 
energy  and  projected  to  an  idealized,  near-ground  image 
plane.  An  atmospheric  model  samples  the  near-ground 
image  and  simulates  transmission  effects  to  the  height  of 
the  sensor.  A sensor  model  simulates  the  effects  of 
optical  blur  and  sensor  noise.  The  resulting  synthetic 


sensor  image  may  be  used  as  if  collected  from  the  field  to 
test  and  improve  automated  target  recognition  algorithms. 

3.  Mesh  Generation  Method 

The  tetrahedral  meshes  used  for  the  moisture  and 
heat  transport  models  are  generated  using  software  that 
takes  advantage  of  the  packages  Triangle1""131  and 
TetGen  [7, 14,1 51  which  are  both  compiled  as  callable 
libraries.  This  mesh  generation  process  has  been 
simplified  to  rapidly  produce  meshes  for  remote  sensing 
simulations  within  the  Countermine  Computational 
Testbed.  The  meshes  can  be  reproduced  quickly, 
reducing  the  need  to  store  large  numbers  of  large  mesh 
files.  This  easy  production  of  meshes  also  greatly  reduces 
the  overall  time  needed  to  run  simulations. 

The  product  of  the  mesh  generation  method  described 
in  this  work  is  a full,  conforming,  three-dimensional  (3D) 
tetrahedral  mesh  of  the  shallow  subsurface  region  under 
investigation.  There  are  no  overset  grids  or  hybrid 
elements.  Objects  are  completely  integrated  into  the 
mesh  and  are  denoted  with  unique  material  identification 
(ID)  for  each  tetrahedron.  The  plants  are  handled  outside 
of  this  scheme.  The  top  of  the  domain,  which  represents 
the  ground  surface,  has  a realistic  topography  based  on 
Light  Detection  and  Ranging  (LIDAR)  data.  These  data 
come  in  the  form  of  (X,Y,Z)  values  and  come  from  field 
measurements,  although  the  data  could  be  generated  using 
geostatistics  or  by  some  other  suitable  method. 
Occasionally,  the  LIDAR  data  are  incomplete  due  to 
“shadow”  effects  from  plants  or  other  obstacles.  In  this 
case,  the  holes  are  patched  using  information  from  the 
immediate  region  around  the  hole. 

Constraints  on  the  maximum  tetrahedral  volume  and 
maximum  area  of  surface  triangles  are  supplied  by  the 
user  at  runtime. 

The  remainder  of  this  section  describes  in  detail  the 
mesh  generation  steps  in  the  order  of  operation. 

3.1.  Top  Surface  Mesh  Generation 

The  surface  of  the  domain,  Figure  2,  was  generated 
using  Shewchuk’s  Triangle  code  guided  by  the  bounding 
region.  After  the  flat  mesh  was  produced,  the  heights,  Z 
values,  were  assigned  based  either  on  LIDAR  input  data 
coming  from  field  measurements  or  data  generated  based 
on  geostatistical  descriptors.  The  resolution  of  the 
LIDAR  imagery  is  generally  higher  than  required  for  the 
mesh  generation.  The  LIDAR  representation  is  coarsened 
to  match  the  specified  mesh  resolution  by  a simple 
smoothing  process.  A simple  averaging  was  used  to 
determine  relative  heights  of  the  surface  nodes,  as  it 
appeared  to  form  more  realistic  topographies  than  inverse 
distance  weighting  techniques.  The  radius  of  influence 
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used  for  this  assignment  was  3 cm,  which  was  based  upon 
the  amount  of  available  input  data. 

Triangle  is  guided  by  an  upper  constraint  on  the  area 
of  the  triangles  supplied  by  the  user  and  a lower  bound  of 
33  degrees  for  the  minimum  angle  constraint  when  it  is 
producing  the  initial  flat  mesh. 

3.2.  Mesh  Remainder  of  Bounding  Domain 

The  vertical  sides  of  the  domain  and  flat  bottom  are 
also  meshed  with  Triangle,  resulting  in  a closed  domain. 
The  same  constraints  on  Triangle  used  for  the  top  of  the 
domain  are  again  used  for  the  remaining  five  sides. 
However,  there  is  one  additional  constraint  on  the 
meshing  of  these  flat  surfaces.  For  each  surface,  no 
additional  nodes  are  allowed  to  be  placed  on  edges  where 
nodes  exist,  i.e.,  the  skeletal  frame  of  the  domain.  This  is 
needed  so  that  the  top,  sides,  and  bottom  may  be  aligned 
to  form  a closed  bounding  region,  without  hanging  nodes. 

Currently,  only  regions  with  a total  of  six  sides  are 
available.  More  general,  polygonal  regions  are  not 
considered. 

3.3.  Object  Insertion,  Mesh  Repair 

One  of  the  most  important  aspects  of  this  process  is 
the  need  to  include  objects  within  the  domain,  e.g.,  rocks 
and  land  mines.  This  procedure  is  handled  automatically. 
These  objects  will  be  fully  integrated  into  the  mesh  and 
therefore  are  added  prior  to  tetrahedralization.  They  may 
be  arbitrarily  placed  and  either  buried,  flush  with,  or 
protruding  from  the  ground  surface.  They  may  not 
protrude  from  the  sides  or  bottom  of  the  domain.  There  is 
no  hard  limit  on  the  number  of  objects  that  can  be 
included.  They  may  be  imported  as  empty,  faceted  shells 
or  may  be  complete  tetrahedral  meshes  with  multiple 
material  regions  within  the  object.  Objects  imported  as 
shells  are  tetrahedralized  by  TetGen  whereas  full  objects 
are  treated  as  holes.  In  the  latter  case,  the  object 
tetrahedra  are  added  back  into  the  domain  afterward. 

Intersections  between  the  objects  and  the  surface  of 
the  domain  are  detected  using  TetGen.  A simple  repair 
operation  is  performed  in  the  immediate  neighborhood  of 
each  object  when  needed  to  “stitch”  the  object  into  the 
surface  mesh.  The  mesh  of  the  object  is  never  altered  but 
a “high  water  mark”  is  identified.  That  is,  nodes  on  the 
object  are  separated  into  two  groups:  those  completely 
above  the  surface  mesh  and  all  others.  Intersecting 
triangles  from  the  surface  mesh  are  removed  and  the 
“apron”  around  the  object  is  then  remeshed  using 
Triangle.  This  apron  reconnects  the  surface  mesh  with 
the  object  shell  along  this  “high  water  mark”.  In  the 
process  the  surface  mesh  may  be  deformed  slightly  and/or 
nodes  may  be  added.  The  object  apron  created  by  the 


reconnection  process  is  not  always  visually  detectable. 
Sensitivity  analysis  shows  that  element  quality  in  the 
apron  regions  has  no  effect  on  the  result  of  simulations 
nor  degrades  performance  of  the  equation  solvers. 

3.4.  Subsurface  Mesh  Generation 

As  stated  before,  the  result  of  the  previous  steps  is  a 
closed  domain  with  objects  that  is  used  as  input  to 
TetGen.  The  only  guidance  given  to  TetGen,  other  than 
the  bounding  region,  is  a volume  constraint  supplied  by 
the  user,  closed  material  regions  denoting  the  objects 
within  the  domain,  and  a request  not  to  subdivide  any 
input  surface  triangle.  Experience  has  shown  that  using 
the  quality  constraints  available  in  TetGen  will 
dramatically  increase  the  time  needed  for  mesh  generation 
for  this  problem.  Therefore,  the  quality  constraint  is 
disabled  in  TetGen  and  the  smoothing  of  the  resulting 
tetrahedra  is  handled  as  a post-process  step. 

3.5.  Tetrahedral  Mesh  Smoothing 

Because  the  mesh  quality  constraints  are  disabled  in 
TetGen  for  this  work,  there  is  no  implied  or  explicit 
guarantee  on  element  quality.  As  a final  step  then,  the 
output  mesh  is  post-processed  with  a mesh-smoothing 
algorithm  to  ensure  quality  tetrahedra  in  the  final  mesh. 
In  this  smoothing  procedure,  the  internal  nodes  of  the 
domain  may  be  moved;  however,  the  connectivity  of  the 
elements  is  not  altered.  The  smoothing  process  employs 
the  element  quality  metric  and  element  size  metric 
detailed  in  References  1 and  2.  The  shape  metric  has 
been  shown  to  be  equivalent  to  other  tetrahedron  quality 
metrics,  as  it  correctly  identifies  all  types  of  poorly 
shaped  elements,  i.e.,  “sliver”  elements,  “needle” 
elements,  “flat”  elements,  etc.  In  general,  “flat”  elements 
tend  to  appear  in  this  process  more  often  than  all  of  the 
other  types  of  poorly  formed  elements  combined.  There 
were  few  of  these,  and  they  tended  to  be  near  the  top 
surface  of  the  domain.  Linking  the  variable  topography 
of  the  surface  to  the  relatively  regular  sides  and  bottom  of 
the  domain  most  likely  caused  these  malformed  elements 
during  the  Delaunay  mesh  generation  process. 

The  size  metric  quantifies  the  difference  between  the 
size  of  a given  element  and  some  optimum  size.  The 
average  element  volume  is  used  throughout  the  domain  as 
the  optimum  in  this  work. 

A variety  of  optimization  techniques  have  been 
employed  in  the  smoothing  step.  Primarily,  the  steepest 
descent  and  Broyden-Fletcher-Goldfarb-Shanno 
algorithms  were  used  for  the  optimization  problem  and 
were  implemented  as  described  in  Reference  9.  The 
gradients  of  the  objective  function  are  derived  and 
computed  explicitly.  Generally,  a static  number  of 
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iterations  (50)  are  taken,  with  a goal  of  improving  the 
final  mesh  as  opposed  to  meeting  a target  quality 
constraint. 

3.6.  Identifying  Subsurface  Soil  Regions 

Subsoil  regions  can  be  assigned  as  a post-process 
step  when  desired.  A distribution  of  soil  regions  can  be 
derived  using  site-specific  geostatistical  information,  i.e., 
borehole  data,  and  the  program  T-PROGS'3-61.  The 
output  of  T-PROGS  is  a Cartesian  grid  with  material  IDs 
distributed  according  to  variogramgenerated  Markov 
chains.  The  program  HYSSOP1161  creates  smoothed 
boundaries  between  the  subsoil  material  regions  using  a 
string-of-pearls  constraint-based  smoothing  algorithm. 
Each  tetrahedron  in  the  mesh  is  assigned  a material  ID 
using  its  centroid  to  determine  its  encompassing  material 
region. 

3.7.  SceneGen  Modeler 

A scene  development  GUI  (graphical  user  interface) 
has  also  been  developed  in  support  of  the  Countermine 
Computational  Testbed.  The  purpose  of  this  modeler  is  to 
orient  and  visualize  the  input  data  prior  to 
tetrahedralization.  SceneGen  allows  the  user  to  visualize 
the  LIDAR  data,  supplied  bounding  region,  and  object 
shells  at  their  relative  locations.  It  also  allows 
visualization  of  vegetation  and  other  energy  occluding 
objects  used  in  the  ray  caster  and  vegetation  models.  This 
ability  to  view  all  of  the  objects  in  the  scene 
simultaneously  allows  the  user  to  plan  and  verily  the 
scene  carefully.  Other  capabilities  of  the  modeler  include 
the  ability  to  move  and  orient  the  objects  within  the  scene. 
The  modeler  can  then  output  the  newly  edited  scene  in 
formats  for  mesh  generation  and  ray  casting. 

A screen  shot  of  the  SceneGen  modeler  can  be  seen 
in  Figure  5.  In  this  scene  both  plants  and  objects  are 
being  moved  within  the  scene  prior  to  mesh  generation. 

4.  Results 

All  tetrahedral  meshes  shown  in  this  paper  were 
created  using  a single  core  on  one  processor  on  a Cray 
XT3  architecture  with  2 GB  of  memory.  In  addition  to 
the  Cray  XT3  architecture,  this  automated  mesh  generator 
package  has  also  been  built  on  the  SGI  Origin  3900  as 
well  as  desktop  machines  running  Linux  and  Mac  OS  X. 
However,  it  is  currently  executed  most  often  on  the  Cray 
XT3  architecture  because  of  the  speed  and  amount  of 
available  memory  of  this  machine.  It  can  be  run  on  the 
front  or  back  end  of  the  XT3. 

At  present  at  least  seven  users  of  varying  educational 
backgrounds  are  creating  meshes.  Over  the  past  year. 


these  seven  users  have  created  approximately  50  meshes 
that  have  been  used  in  hundreds  of  heat  transport 
simulations.  Because  the  typical  user  needs  to  change 
only  the  number  and  (X,Y,Z)  locations  of  objects,  the 
time  needed  for  training  new  users  is  generally  well  under 
an  hour. 

The  total  time  needed  to  create  the  mesh  in  Figure  1 
was  less  than  one  hour,  most  of  which  was  spent  in  the 
planning  of  object  locations.  There  are  approximately  1.5 
million  tetrahedral  elements  and  270  thousand  nodes. 
The  total  ‘computational  time  was  nine  minutes  with  the 
majority  of  that  time  spent  in  the  mesh,  smoothing  out  the 
step  (approximately  3-4  minutes).  In  this  manner,  the 
objects  could  be  moved  (or  others  could  be  added)  and  the 
new  mesh  could  be  generated  in  well  under  an  hour. 
Therefore,  the  testing  of  new  scenarios  with  a similar 
configuration  can  be  achieved  with  minimal  effort. 

The  surface  of  the  domain  is  shown  in  Figure  2.  This 
surface  had  57  thousand  triangles  prior  to  the  insertion  of 
objects,  and  the  resulting  topography  can  clearly  be  seen 
in  the  figure.  There  are  a total  of  four  rocks  in  the  final 
domain,  which  can  be  seen  in  the  full  mesh  in  Figure  1. 
The  colors  are  determined  by  material  type.  A closer 
view  of  two  of  the  rocks  can  be  seen  in  Figure  3.  The 
object  intersection  routines  have  seamlessly  integrated  the 
rocks  into  the  scene. 

A separate  mesh  was  created  using  roughly  the  same 
area,  and  the  result  of  a full  simulation  can  be  seen  in 
Figure  4.  In  this  figure  the  colors  correspond  to 
temperatures  and  there  is  preferential  loading  in  the 
direction  of  the  sun,  with  shadows  being  cast  by  the 
objects.  This  demonstrates  the  interactions  of  the  solar 
model  with  the  vegetation  model  and  the  heat  transport 
model  in  the  shallow  subsurface. 

5.  Future  Work 

Large  field  scale  (1  km  x 1 km)  simulations  are  being 
planned  for  the  near  future.  The  meshing  approach 
described  in  this  paper  is  near  the  memory  limitations  of  a 
single  processor  on  the  Cray  XT3  being  used.  Initially, 
larger  domains  will  be  created  in  a patchwork  fashion 
with  matching  vertical  boundaries  between  the 
subdomains.  Simple  scalability  will  be  achieved  by  using 
a separate  processor  for  each  subdomain. 

Currently,  tetrahedral  element  volume  is  determined 
by  the  resolution  needed  on  the  surface  of  the  domain. 
However,  the  resulting  number  of  tetrahedra  is  quite 
large.  The  number  of  elements  produced  could  be 
reduced  by  grading  the  mesh  in  the  negative  z-direction, 
especially  below  the  depth  of  diurnal  temperature 
fluctuations. 
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6.  Conclusion 

The  process  of  rapid,  unstructured  tetrahedral  mesh 
generation  for  remote  sensing  simulations  has  been 
achieved  on  a variety  of  computational  architectures.  The 
final  mesh  can  include  multiple  objects,  each  at  an 
arbitrary  location  within  the  domain  and  the  mesh 
generally  can  be  reproduced  in  much  less  than  an  hour. 
Individual  scenes  can  be  planned  in  detail,  allowing 
objects  to  be  carefully  placed,  thereby  greatly  increasing 
the  number  of  scenarios  that  can  be  tested. 
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Figure  1.  The  full  3D  mesh  with  272k  nodes  and  1.5M 
tetrahedra  shown  prior  to  the  addition  of  plants 


Figure  2.  Surface  mesh,  prior  to  object  insertion,  with 
57,201  triangles  covering  a square  domain,  10  meters 
on  a side 
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Figure  3.  Close-up  of  two  rocks.  Notice  that  the 
aprons  around  the  rocks  are  blending  into  the  full 
mesh. 


Figure  4.  Region  with  many  objects,  two  types  of 
plants  in  various  locations,  and  the  resulting 
temperatures  during  the  course  of  a simulation 


Figure  5.  Screen  Shot  of  SceneGen  Modeler  Software. 

In  this  figure,  there  are  several  objects  as  well  as 
several  plants  being  incorporated  into  the  scene.  The 
highlighted  objects  are  in  the  process  of  being  moved. 


Figure  6.  A view  along  the  surface  of  a mesh  with 
multiple  plants  and  objects.  Notice  the  rock  in  the 
scene  which  has  been  blended  into  the  surface  mesh. 
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